make_scale_coefs <- function(.df, .dep_var, ...){
  
  #dep vars scaling
  .dep_var = enquo(.dep_var)
  .df =
    .df %>%
    mutate(outcome = !!.dep_var) %>%
    mutate(scaled_dv_outcome = scale(outcome),
           scaled_iv_pop = scale(log10(pop)))
  
  #independent vars scaling
  .ind_vars = enquos(...)
  for(i in .ind_vars){
    .df =
      .df %>%
      mutate(ind_var = scale(!!i))
    
    names(.df)[names(.df) == "ind_var"] = paste0("scaled_iv_", quo_name(i))
    
  }
  
  
  form_reg <-
    as.formula(
      paste0("scaled_dv_outcome ~ scaled_iv_pop +",
             tidyselect::vars_select(names(.df), starts_with('scaled_iv_', ignore.case = TRUE)) %>%
               paste0(., collapse = " + ")))
  
  lm(formula = form_reg, data = .df) %>%
    tidy() %>%
    return()
  
}

brazil_metro <- 
  brazil_data %>%
  filter(urb == 1, metro_area == 1)

brazil_non_metro <- 
  brazil_data %>%
  filter(urb == 1, metro_area == 0)



doc_mod_brazil_met <- make_scale_coefs(.df = brazil_metro, .dep_var = prop_doctos, ... = pop_growth, density, div_index, avg_salary_mo_2010, ill_rate_2010, elec_index, align) %>%
  filter(term == "scaled_iv_pop") %>% mutate(country = "Brazil") %>% mutate(outcome = "docto") %>% mutate(type = "Metro")
teacher_mod_brazil_met <- make_scale_coefs(.df = brazil_metro, .dep_var = teachers_prop, ... = pop_growth, density, div_index, avg_salary_mo_2010, ill_rate_2010, elec_index, align) %>%
  filter(term == "scaled_iv_pop") %>% mutate(country = "Brazil") %>% mutate(outcome = "teacher") %>% mutate(type = "Metro")
comm_mod_brazil_met <- make_scale_coefs(.df = brazil_metro, .dep_var = prop_comm_health, ... = pop_growth, density, div_index, avg_salary_mo_2010, ill_rate_2010, elec_index, align) %>%
  filter(term == "scaled_iv_pop") %>% mutate(country = "Brazil") %>% mutate(outcome = "health") %>% mutate(type = "Metro")
schools_mod_brazil_met <- make_scale_coefs(.df = brazil_metro, .dep_var = schools_prop, ... = pop_growth, density, div_index, avg_salary_mo_2010, ill_rate_2010, elec_index, align) %>%
  filter(term == "scaled_iv_pop") %>% mutate(country = "Brazil") %>% mutate(outcome = "schools") %>% mutate(type = "Metro")
elec_mod_brazil_met <- make_scale_coefs(.df = brazil_metro, .dep_var = electr_share_2010, ... = pop_growth, density, div_index, avg_salary_mo_2010, ill_rate_2010, elec_index, align) %>%
  filter(term == "scaled_iv_pop") %>% mutate(country = "Brazil") %>% mutate(outcome = "elect") %>% mutate(type = "Metro")
water_mod_brazil_met <- make_scale_coefs(.df = brazil_metro, .dep_var = piped_water_share_2010, ... = pop_growth, density, div_index, avg_salary_mo_2010, ill_rate_2010, elec_index, align) %>%
  filter(term == "scaled_iv_pop") %>% mutate(country = "Brazil") %>% mutate(outcome = "water") %>% mutate(type = "Metro")

doc_mod_brazil_nonmet <- make_scale_coefs(.df = brazil_non_metro, .dep_var = prop_doctos, ... = pop_growth, density, div_index, avg_salary_mo_2010, ill_rate_2010, elec_index, align) %>%
  filter(term == "scaled_iv_pop") %>% mutate(country = "Brazil") %>% mutate(outcome = "docto") %>% mutate(type = "non-Metro")
teacher_mod_brazil_nonmet <- make_scale_coefs(.df = brazil_non_metro, .dep_var = teachers_prop, ... = pop_growth, density, div_index, avg_salary_mo_2010, ill_rate_2010, elec_index, align) %>%
  filter(term == "scaled_iv_pop") %>% mutate(country = "Brazil") %>% mutate(outcome = "teacher") %>% mutate(type = "non-Metro")
comm_mod_brazil_nonmet <- make_scale_coefs(.df = brazil_non_metro, .dep_var = prop_comm_health, ... = pop_growth, density, div_index, avg_salary_mo_2010, ill_rate_2010, elec_index, align) %>%
  filter(term == "scaled_iv_pop") %>% mutate(country = "Brazil") %>% mutate(outcome = "health") %>% mutate(type = "non-Metro")
schools_mod_brazil_nonmet <- make_scale_coefs(.df = brazil_non_metro, .dep_var = schools_prop, ... = pop_growth, density, div_index, avg_salary_mo_2010, ill_rate_2010, elec_index, align) %>%
  filter(term == "scaled_iv_pop") %>% mutate(country = "Brazil") %>% mutate(outcome = "schools") %>% mutate(type = "non-Metro")
elec_mod_brazil_nonmet <- make_scale_coefs(.df = brazil_non_metro, .dep_var = electr_share_2010, ... = pop_growth, density, div_index, avg_salary_mo_2010, ill_rate_2010, elec_index, align) %>%
  filter(term == "scaled_iv_pop") %>% mutate(country = "Brazil") %>% mutate(outcome = "elect") %>% mutate(type = "non-Metro")
water_mod_brazil_nonmet <- make_scale_coefs(.df = brazil_non_metro, .dep_var = piped_water_share_2010, ... = pop_growth, density, div_index, avg_salary_mo_2010, ill_rate_2010, elec_index, align) %>%
  filter(term == "scaled_iv_pop") %>% mutate(country = "Brazil") %>% mutate(outcome = "water") %>% mutate(type = "non-Metro")

br_metro <-
  bind_rows(comm_mod_brazil_nonmet, schools_mod_brazil_nonmet, doc_mod_brazil_nonmet, teacher_mod_brazil_nonmet, elec_mod_brazil_nonmet, water_mod_brazil_nonmet,
          comm_mod_brazil_met, schools_mod_brazil_met,doc_mod_brazil_met, teacher_mod_brazil_met, elec_mod_brazil_met, water_mod_brazil_met) %>%
  mutate(label = case_when(str_detect(outcome, "schools") ~ "Schools (per 1000)",
                           str_detect(outcome, "teacher") ~ "Teachers (per 1000)",
                           str_detect(outcome, "water") ~ "Piped Water (%)",
                           str_detect(outcome, "elect") ~ "Electricity Access (%)",
                           str_detect(outcome, "docto") ~ "Doctors (per 1000)",
                           str_detect(outcome, "health") ~ "Health Centers (per 1000)",
                           TRUE ~ NA_character_)) %>%
  mutate(rank = case_when(str_detect(outcome, "schools") ~ "2",
                          str_detect(outcome, "water") ~ "9",
                          str_detect(outcome, "elect") ~ "8",
                          str_detect(outcome, "teacher") ~ "6",
                          str_detect(outcome, "docto") ~ "5",
                          str_detect(outcome, "health") ~ "3",
                          TRUE ~ NA_character_)) %>%
  dplyr::select(estimate, std.error, label, type, rank) %>%
  add_row(estimate = NA_real_, std.error = NA_real_, label = "Divisible Public Goods:                             ", type = "Metro", rank = "1") %>%
  add_row(estimate = NA_real_, std.error = NA_real_, label = "Personnel:                                          ", type = "Metro", rank = "4") %>%
  add_row(estimate = NA_real_, std.error = NA_real_, label = "Indivisible Public Goods:                           ", type = "Metro", rank = "7") %>%
  mutate(rank = as.numeric(rank)) %>%
  ggplot(aes(x=estimate, y = reorder(label, -rank), color = type)) +
  geom_point(position = position_dodge(width = 1)) +
  geom_errorbarh(aes(xmin=estimate-1.96*std.error, xmax = estimate+1.96*std.error), position = position_dodge(width = 1), height = 0.1) +
  geom_vline(xintercept = 0, linetype = "dotted") +
  #facet_wrap(country~.) +
  theme_bw() +
  scale_colour_grey() +
  theme(panel.grid.minor = element_blank(), 
        panel.grid.major.x = element_blank(),
        axis.line.y.left = element_blank(),
        axis.ticks.y = element_blank(),
        legend.position = "bottom",
        strip.background = element_blank(),
        legend.title = element_blank(),
        axis.line = element_line(colour = "black"),
        panel.border = element_blank(),
        axis.title.y = element_blank()) +
  xlab("Standardized Coefficient")

ggsave("./_4_outputs/figure_oa_iii.pdf", plot = br_metro, width = 5, height = 3)
